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Abstract 

Given n independent random variables Xi,X2, ■ ■■,X n and an integer C, we study the fundamental 
problem of computing the probability that the sum X — Xi + X2 + ... + X n is at most C. We assume 
that each random variable Xi is implicitly given by an oracle which, given a input value k, returns 
the probability Xi < k. We give the first deterministic fully polynomial-time approximation scheme 
(FPTAS) to estimate the probability up to a relative error of 1 ± e. Our algorithm is based on the idea 
developed for approximately counting knapsack solutions in [Gopalan et al. FOCS11]. 

1 Introduction 

Wc study the following fundamental problem: The input consists of n independent (not necessarily identically 
distributed) random variables X±, . . . , X n and an integer C. Let their sum be X = X\ + X 2 + ■ . . + X n . 
Our task is to compute the following probability value 

k 



F(C7)=Pr[^A 4 <C7 



Assumptions: We assume that all random variables are discrete and the support of Xi, denoted as supp^, 
consist of only integers. We also assume that each random variable Xi is implicitly given by an oracle 
(denoted as OC) which, given a input value k, returns the probability Xi < k. Note that this assumption is 
weaker than assuming the explicit representation of Xi (listing the probability mass at every point), since we 
need to use binary search with IsuppJ log supp.J calls to the oracle to construct the explicit representation 
of*,. 

It is well known that computing F(C) is #P-hard (see e.g., [9]). The hardness of computing F(C) has 
an essential impact in the area of stochastic optimization as many problems generalize and/or utilize this 
basic problem in one way or another, thus inheriting the #P-hardness. Although we can sometimes use for 
example the linearity of expectation to bypass the difficulty of computing F(C), more than often no such 
simple trick is applicable, especially in the context of risk-aware stochastic optimization where people usually 
pay more attention to the tail probability than the expectation. 

Despite the importance of the problem, surprisingly, no approximation algorithm with provable mul- 
tiplicative factor is known. We note that we can easily obtain an additive PTAS for this problem via the 
Monte-Carlo method: One generates K independent samples X"\ i = 1,2, K of X and use the empirical 
sum 

i=l 
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1 So, the running time of a polynomial time algorithm should be a polynomial of log C. 
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where /(•) is the indicator function. By standard Chernoff bound, one can easily see that with K = poly(l/e) 
samples, the estimation is within an additive error e from the true value with a constant probability (see 
e -g-i j!4j). To get a reasonable multiplicative approximation factor, we need to set the value of e at the order 
of F(C), so the number of samples needs to be poly(l/F(C)), which can be exponentially large, when F(C) 
is exponentially small 0. 

In this paper, we propose a fully polynomial approximation scheme (FPTAS) for counting F(C). For 
ease of notation, we use (1 ± e)F(C) to denote the interval [(1 — e)F(C), (1 + e)F(C)]. Recall that we say 
there is an FPTAS for the problem, if for any fixed constant e > 0, the algorithm can produce a value F 
such that F € (1 ± e)F(C) in polynomial (w.r.t. the length of the binary encoding of the input) time (See 
e.g., El). 

Theorem 1.1. We are given n independent random variables X\, . . . , X n distributed on support [N], a 
positive integer C, a fixed positive constant e (e > 0). Suppose there is an oracle for each Xi, which, 
upon two input integers (ni,ri2), returns the value Pr[ni < Xi < n?\. There is an FPTAS for estimating 
F >r E"=i^ ^ C] an d th e running time is 0(p- log(i) 2 log(niV)) where q = Yii mhizesuppCXi) Pr[Xj = a?]. 

In the above theorem, we assume the binary coding of the value Pr[A"i = x] is polynomially bounded 
for any x £ supp(Ai). Hence, logg is polynomial. Note that in the theorem, we can assume TV < C + 1 by 
moving all probability mass exceeding C + 1 to point C + 1 without affecting the final answer. 

1.1 Related Works 

There is a large body of work on estimating or upper/lower-bounding the distribution of the sum of indepen- 
dent random variables. See e.g., [TJ [18J [12l El [13] . Those works are based on analytic numerical methods (e.g., 
Edgeworth expansion, saddle point method) which either require specific families of distributions and/or do 
not provide any provable multiplicative approximation guarantees. 

Our problem is a generalization of the counting knapsack problem. For the counting knapsack problem, 
Morris and Sinclair [15] obtained the first FPRAS (fully-polynomial randomized approximation scheme). 
Dyer [4] provided a completely different FPRAS based on dynamic programming. The first deterministic 
FPTAS is obtained by Gopalan et al.[7]. 

Our problem is also closely related to the threshold probability maximization problem (see a general 
formulation in |llj). In this problem, we are given a ground set of items. Each feasible solution to the 
problem is a subset of the elements satisfying some property (this includes problems such as shortest path, 
minimum spanning tree, and minimum weight matching). Each element b is associated with a random 
weight Xb. Our goal is to to find a feasible set S such that Pr Ebes^>] i s maximized. There is a large 
body of literature on the threshold probability maximization problem, especially for specific combinatorial 
problems and/or special distributions. For example, Nikolova, Kelner, Brand and Mitzenmacher [17] studied 
the corresponding shortest path version for Gaussian, Poisson and exponential distributions. Nikolova [16] 
extended this result to an FPTAS for any problem with Gaussian distributions, if the deterministic version of 
the problem has a polynomial time algorithm. The minimum spanning tree version with Gaussian distributed 
edges has also been studied in [5]. For general discrete distributions, Li and Deshpande [TU] obtained an 
additive PTAS if the deterministic version of the problem can be solve exactly in pseudopolynomial time. 
Very recently, Li and Yuan jllj further generalized this result to the class of problems for which the multi- 
objective deterministic version admits a PTAS. 

Our problem is also closely related to the fixed set version of the stochastic knapsack problem. In this 
problem, we are given a knapsack of capacity C and a set of items with random sizes and profits. Their goal 

2 In certain application domains such as risk analysis, small probabilities (often associated with catastrophic losses) can be 
very important. 
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is to find a set of items with maximum total profit subject to the constraint that the overflow probability is 
at most a given parameter 7. Kleinberg, Rabani and Tardos [9] first considered the problem with Bernoulli- 
type distributions and provided a polynomial-time 0(log l/7)-approximation. Better results are known for 
specific distributions, such as exponentially distributions jBJ, Gaussian distributions [5J [T^]. For general 
discrete distributions, bi-criteria additive PTASesE) arc known via different techniques [2l [TU1 [TTj . 

2 Algorithm 

Our FPTAS is based on dynamic programming. In Section 12. 1[ we provide the recursion of the dynamic 
program, which is largely based on the idea developed in [7], with some necessary adaptations. However, 
since the support of each random variable can be exponentially large, it is not immediate clear how the 
recursion can be implemented efficiently given the oracles. In Section [2~2l we address this issue. 

2.1 The Dynamic Program 

We first notice that PrE*. =1 Xj < C] is a nondecreasing function of C. We consider an inverse function 
r(i,a) : {0,1, ...,n} x R — > R U {±00}, which is defined to be the smallest C such that the probability 
PrEj"=i Xj < C] > a. It is easy to sec that r(i, a) is nondecreasing in a. We need the following simple 
lemma. We omit the proof, which is straightforward. 

Lemma 2.1. 

n 

Pr[^ X t <C} = max{a : r(n, a) < C}. 
3=1 

The following recursion is very important to us. 

Lemma 2.2. Suppose di S supp(Wi), p(di) = Pr[X^ = di], 7 : supp(X;) — > [0,1]. We have the following 
recurrence 

r(i, a) = minmax{r(i — 1, 7(<i)) + di} (2) 

7 di 

subject to 

Y l{di)p{di) > a (3) 

di Esupp( Xi ) 

Intuitively, suppose that j(di) represents the probability Pr[X)}=i Xj < C — di] for some C. <j3j) would 
enforce that Pr[^}=i Xj < C] = X)diesup P ( xA P(di)j{di) > a. The formal proof is as follows. 

Proof. We first prove that for every choice of 7 : supp(X,) — s> [0, 1] (a list of 7(2^) for xi € supp(Xj)), the 
quantity C ~ rnaxd fcgsU pp(x i ){'K* — lj7(^fc)) + ^fe} > r(i,a). Due to the independence, We can see that 

i i—1 

Pr[]T Xj <C']= Pr E X i ^ C '- d *l Pr [*< = d *l 

3=1 d & esupp(Xi) j=l 

By the definition of C , we know that C — dk > r(i — l,j(dk)) for any dk £ supp(Wi), thus Pr[^*Cj Xj < 
C'-d k ] > 7 (4). So 

i 

Pv[Y,Xj<C'}> ]T l(d k )p{d k )>a 

j=l d k esupp(Xi) 
3 The overflow probability constraint may be violated by an additive factor e for any constant e > 0. 
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Algorithm 1 The Dynamic Program for Computing Pr[^ X t < C] 

1: Set T[l,j] = t(1,Q-J) for all j > 0. 

2: Set Q = 1 + i^i±£i and s = r i ogQ I] . 

3: for i = 2 — !> k do 

4: for j : = — > s do 

5: Set T(i,j) = min A max d , T(i - 1, j + [log Q {^)\) + d, 

6: end for 

7: end for 

8: Let f = min{j : T[nJ] < C}. Output F' := Q~i' +1 . 



Second, we prove that there exists a choice of 7' such that C = maxd k {r(i — l,7'(dfc)) + dk} < r(i,a). 

Let 

i-l 

y^HPr^X,- <r(i,a)-4] ( 4 ) 
3=1 

It is easy to see that 7' satisfies By definition, r(i — I, 7 / (d/ ; )) < r(i, a) — dj. for all dj. € supp(Xj), which 

leads to max dfc {r(i — l,7'(dfc)) + d k } < r(i,a). This completes the proof. □ 

For ease of notation, we let A(dj)et = 7(d i )p(d i ). The recursion can be written as: 

r(i, a) = minmax{r(i — 1, a) + di} subject to Mdi) > 1 (5) 

A di v(di) ^-^ 

V ' dj6supp(Xi) 

It is not clear how recursion ([5]) can be efficiently executed since the second argument a is a continuous 
variable. We use the following way to discretize it: Let Q = 1 + ln ^^ . Note that Q is slightly larger than 
1. We use q to denote I7 ie r n i m i n d i esupp(jf i ) Pr[-Xi = di], which is clearly a lower bound of the answer. Let 

a = riog Q 11 = o(aiog 2 1). 

We define recursively the function T : {0, 1, n} x {0, 1, s}->lU {±00} as follows: 

T(i,j)=rrunmaxT(i-l,j+Llog Q (^|)j)+d i subject to ^ \{di) > 1 (6) 

di£supp(Xi) 

Comparing ([5]) and ([S]), the similarity suggests that T(i,j) is the approximate version of r(i, Q~ 3 ). The next 
lemma formalizes this idea. 

Lemma 2.3. For all i £ {0, 1, n} and j S {0, 1, s}, we have 

T(i,Q-j)<T(i,j)<r(i,Q-U-V) (7) 

Proof. We prove this by induction. The base case is trivial by the definition of T(0,j). Now assume for the 
statement is true for i — 1 and all j 6 {0, 1, 2, s}. We prove it is also true for i and all j G {0, 1, 2, s}. 
By induction hypothesis, we have that 

T(i - 1, [j + log Q fjgj) < r(i - 1, Q"(^o gQ < r( . _ x Wg-y-oj and 

T(i - 1, b + logQf^J) > r(i - 1, Q-^mi) > r(i - I, ^Q"'). 
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Algorithm 2 Efficient Implementation of Recursion © 

l: Set Left = 0, Right = kN. 

2: while Right > Left do 

3: Set V = [(Left + Right) /2J . 

4: for ni = 1 — y s do 

5: Let p(j,m) = Pr[Xi G (T" — T(i - l,m- 1),T' -T(i- l,m)] (via oracle O,). 

6: end for 

7: if EL=iQ m "Mj,™)>lthen 

8: Set Right = T> 

9: else 

10: Set Left = T" + 1. 

11: end if 

12: end while 

13: Set T(i,j) = T. 



Taking the maximum over di and the minimum over A do not change the direction of the inequalities. 
Combining the above inequalities with ([5]), we complete the proof. □ 

With this lemma, we can approximate the recursion [5] by solving the rccursion[6] The pseudo-code of the 
dynamic program is provided in Algorithm 1. It is not hard to show the output of the algorithm is a good 
approximation of the true probability. 

Lemma 2.4. The output F' of the algorithm is a (1 ± e) -approximation of F(C) = Pr[J^" Xi < C). 

Proof. From the choice of j', we know that T(n,j') < C < T(n,j' — 1). According to Lemma \2. 31 we have 
T(n,Q- j ') < C < r(?7,Q"-j" +1 ). Consequently F = Pr[£" =1 JTj < C] G [Q~ j ' , Q n ~ j ' +1 ]. By outputting 
F' = Q~ j ' +1 , the approximation ratio can bounded as 1 - e < Q -(n+1 ) < < Q n+1 = 1 + e. □ 

2.2 An Efficient Implementation using Binary Search 

In this subsection, we show how to implement the recursion ([6]) in polynomial time. Suppose we have 
already computed T(i',j') for all i' < i and < j' < s and are computing the value T(i,j). Our approach 
is based on a binary search on the range of T(i,j). For every guess X", we can decide whether T(i,j) < T' 
efficiently by using the criterion in Lemma 12.51 The pseudocode is provided in Algorithm 2. 

Lemma 2.5. We define X'(d,) = max{A(d i ) | T(i - 1, j + L lo gQ(l^y)J ) + d i < T '}- Then > T (hj) < T ' */ 
and only i/Ed ie supp(X0 ~ L 

Proof. Suppose Sd esupp(X ) ^'(^0 — ^- This means there exists a choice of A such that 

maxT(i — Llog Q (^|)j) + d, < T', 

and therefore T(i,j) < T' by ©. On the other hand, suppose X^esuppfx-) A'(dj) < 1. Clearly, any other 
choice of A(di) should be smaller than X'(di) in order to obtain T(i,j) < T', which is not possible. □ 

From Lemma T2.51 we can see that if the |supp(Ai)| is bounded by a polynomial, we can decide whether 
T(i,j) < T' in polynomial time. However, |supp(Xi)| can be exponentially large and only an oracle Oi is 
provided. We resolve this issue by noticing that the support of each random variable can be divided into s 
segments such that the sum of X'(di) values over each segment can be computed efficiently. In particular, 
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we divide [N] into segments (V -oo,T'- T(i - 1,0)] n [TV], (V - T(i - 1,0), V - T(i - 1, 1)] n [TV], {T - 
T(i — 1, s — 1), T' — T(i — 1, s)] n [JV]. In each segment, A'(dj) can be computed according to the following 
lemma. 

Lemma 2.6. Suppose di £ (T' — T(i — l,m — 1),T" — r(? — l,m)] for some m £ [l,s]. Then 

V(di) = Q m - j P (di) 

Proof. For any £ (T' - T(i - 1, m - 1), T' - T(t - 1, m)], by definition of A', we have that 

So A'(di) = Q m ~ip(di) and we complete the proof. □ 

As a result, it is sufficient to ask the Oi for the probability value Pr[X; £ (T" — T(i — 1, m — 1), T'— T(i — 
l,m)]] and we can compute Z)^>T'-T(i-i,m-i) ^'(^) sim ply by 

T'-T(i-l,m) T'-T(i-l,m) 

£ A'(di)= ^ Q m - J *Pr[X i e(T'-T(i-l,m-l) j r'-r(i-l,m)]]. 

di>T'-T(i-l,m-l) di>T'-T(i-l,m-l) 

Lemma 2.7. Suppose each query to the oracle takes O(l) time. The optimization with respect to X,di in 
equation |0J) runs in time 0((^) 2 (log log 2 (niV)) 

Proof. The binary search takes log 2 (nN) iterations, while in each iteration we compute A' for at most s 
segments. For each segment, we query an oracle Oi, which takes one unit of time, to obtain the probability 
of the random variable Xi being in the segment. Q m_3 can be computed in advance and its contribu- 
tion to running time is neglected. So the implementation of each recursion step runs in 0(slog 2 (nN) = 
0(2 1og(i)log(niV)). □ 

Since the number of goal states in the dynamic program is 0(ns), according to Lemma 12.71 the total 
running time is 0(^r l°g(^) 2 log(n./V)). This completes the proof of Theorem ll.il 
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